# Adjust tex table spacing
adjust_spacing <- function(tex_code){
  tex_code <- gsub("\\\\begin\\{tabular\\}","\\\\renewcommand{\\\\arraystretch}{0.75}\n\\\\begin{tabular}",tex_code)
  return(tex_code)
}

# Thinning function
pdrop <- function(x, lower = 0.4, upper = 0.6, na.rm = TRUE) {
  qs <- quantile(x, probs = c(lower, upper), na.rm = na.rm)
  x[x > qs[1] & x < qs[2]]
}

# Extract estimates from SAR/CAR model output
model_scar <- function(y,x,logy=TRUE){
  if(logy==TRUE){fhead<-paste0("log(",y,"+1)")}else{fhead<-y}
  out <- as.formula(paste0(fhead,"~log(",x,"+1) + TROIKA_DIST  + splines::bs(LONG)*splines::bs(LAT)+as.factor(RAY_NAME)+as.factor(BLOK_ZONE)"))
  return(out)
}
scar_coef <- function(modz,tvar,fez=NULL){
  coef_t <- summary(modz)$coefficients$fixed%>%.[grep(tvar,rownames(.)),] 
  pr2 <- summary(modz)$pseudoR2
  nz <- modz$n
  rmsez <- sqrt(mean((modz["residuals"][[1]]["response"][[1]])^2))
  aicz <- stats::AIC(modz)
  sp_covz <- summary(modz)$coefficients$spcov
  ynamz <- attr(modz["terms"][[1]],"dataClasses")%>%names()%>%.[1]
  xnamz <- attr(modz["terms"][[1]],"dataClasses")%>%names()%>%.[2]
  if(length(fez)>0){fe1<-fez[1]}else{fe1<-""}
  if(length(fez)>1){fe2<-fez[2]}else{fe2<-""}
  out <- data.table::data.table(
    model=modz$call$spcov_type,
    y=ynamz,
    x=xnamz,
    beta_est=coef_t[1,1]%>%as.numeric(),
    beta_se=coef_t[1,2]%>%as.numeric(),
    beta_z=coef_t[1,3]%>%as.numeric(),
    beta_p=coef_t[1,4]%>%as.numeric(),
    fe1=fe1,
    fe2=fe2,
    sp_de=sp_covz["de"]%>%as.numeric(),
    sp_ie=sp_covz["ie"]%>%as.numeric(),
    sp_range=sp_covz["range"]%>%as.numeric(),
    sp_extra=sp_covz["extra"]%>%as.numeric(),
    pseudo_r2=pr2%>%as.numeric(),
    rmse=rmsez,
    aic=aicz%>%as.numeric()%>%round(1),
    n=nz%>%as.integer()
    )
  return(out)
}

# Extract estimates from FRDD
model_rddz <- function(y,x,logy=TRUE){
  if(logy==TRUE){fhead<-paste0("log(",y,"+1)")}else{fhead<-y}
  suppressMessages({
    out <- Formula::Formula(paste0("c(",paste0(fhead,collapse=", "),") ~ splines::bs(d_rayon, 3)   |  0 | log(",x,"+1) ~ I(d_rayon > 0) ")%>%as.formula())
  })
  return(out)
}

# Extract estimates from feols
feols_coef <- function(modz,tvar,fez=NULL,ix=NULL){
  if(!class(modz)%in%"fixest_multi"){modz <- list(modz)}
  out <- lapply(1:length(modz),function(m0){#print(m0)
    coef_t_all <- modz[[m0]]$coeftable 
    coef_t <- modz[[m0]]$coeftable %>%.[grep(tvar,rownames(.)),]
    if(any(class(coef_t)!="numeric")){coef_t<-coef_t%>% .[1,] }
    conv_status <- ifelse(is.null(modz[[m0]]$convStatus),TRUE,modz[[m0]]$convStatus)
    off_status <- !is.null(modz[[1]]$call$offset)
    if("iv"%in%names(modz[[m0]])){
      iv1_f <- fixest::fitstat(modz[[m0]], "ivf")[[1]]["stat"][[1]]
      iv1_fp <- fixest::fitstat(modz[[m0]], "ivf")[[1]]["p"][[1]]
      iv1_wu <- fixest::fitstat(modz[[m0]], "wh")[[1]]["stat"][[1]]
      iv1_wup <- fixest::fitstat(modz[[m0]], "wh")[[1]]["p"][[1]]
    } else {
      iv1_f <- NA_real_
      iv1_fp <- NA_real_
      iv1_wu <- NA_real_
      iv1_wup <- NA_real_
    }
    rmsez <- fixest::fitstat(modz[[m0]],"rmse")[[1]]
    nz <- fixest::fitstat(modz[[m0]],"n")[[1]]
    aicz <- fixest::fitstat(modz[[m0]],"aic")[[1]]
    ar2 <- fixest::fitstat(modz[[m0]],"ar2")[[1]][[1]]
    pr2 <- fixest::fitstat(modz[[m0]],"pr2")[[1]][[1]]
    if(!is.na(ar2)){pr2 <- NA}
    ynamz <- gsub("lhs\\: ","",names(modz)[[m0]])%>%stringr::str_split("\\;",simplify=TRUE)%>%.[1,1]
    xnamz <- modz[[m0]]$coeftable%>%rownames()%>%grep(tvar,.,value=TRUE)
    if(!is.null(ix)){
      ixnamz <- modz[[m0]]$coeftable%>%rownames()%>%grep(ix,.,value=TRUE)
      if(length(ixnamz)==1){
        beta_ixl_est <- NA_real_
        beta_ixl_se <- NA_real_
        beta_ixl_z <- NA_real_
        beta_ixl_p <- NA_real_
        beta_ixu_est <- coef_t_all[ixnamz,][1]%>%as.numeric()
        beta_ixu_se <- coef_t_all[ixnamz,][2]%>%as.numeric()
        beta_ixu_z <- coef_t_all[ixnamz,][3]%>%as.numeric()
        beta_ixu_p <- coef_t_all[ixnamz,][4]%>%as.numeric()
      }
      if(length(ixnamz)>1){
        beta_ixl_est <- coef_t_all[ixnamz,][1,][1]%>%as.numeric()
        beta_ixl_se <- coef_t_all[ixnamz,][1,][2]%>%as.numeric()
        beta_ixl_z <- coef_t_all[ixnamz,][1,][3]%>%as.numeric()
        beta_ixl_p <- coef_t_all[ixnamz,][1,][4]%>%as.numeric()
        beta_ixu_est <- coef_t_all[ixnamz,][2,][1]%>%as.numeric()
        beta_ixu_se <- coef_t_all[ixnamz,][2,][2]%>%as.numeric()
        beta_ixu_z <- coef_t_all[ixnamz,][2,][3]%>%as.numeric()
        beta_ixu_p <- coef_t_all[ixnamz,][2,][4]%>%as.numeric()
      }
    } else {
      beta_ixl_est <- beta_ixl_se <- beta_ixl_z <- beta_ixl_p <- beta_ixu_est <- beta_ixu_se <- beta_ixu_z <- beta_ixu_p <- NA_real_
    }
    if(length(fez)>0){
      fe_ <- fez %>% as.list() %>% data.table::as.data.table() %>% data.table::setnames(paste0("fe",1:length(fez)))
      out_ <- data.table::data.table(
              model=modz[[m0]]["method"][[1]],
              y=ynamz,
              x=xnamz[1],
              ix=ifelse(is.null(ix),NA_character_,ix),
              beta_est=coef_t[1]%>%as.numeric(),
              beta_se=coef_t[2]%>%as.numeric(),
              beta_z=coef_t[3]%>%as.numeric(),
              beta_p=coef_t[4]%>%as.numeric(),
              beta_ixl_est = beta_ixl_est,
              beta_ixl_se = beta_ixl_se,
              beta_ixl_z = beta_ixl_z,
              beta_ixl_p = beta_ixl_p,
              beta_ixu_est = beta_ixu_est,
              beta_ixu_se = beta_ixu_se,
              beta_ixu_z = beta_ixu_z,
              beta_ixu_p = beta_ixu_p,
              converge=conv_status,
              offset=off_status,
              iv1_f=iv1_f,
              iv1_fp=iv1_fp,
              iv1_wu=iv1_wu,
              iv1_wup=iv1_wup) %>% dplyr::bind_cols(fe_) %>% dplyr::bind_cols(
                data.table::data.table(
                  adj_r2=ar2%>%as.numeric(),
                  pseudo_r2=pr2%>%as.numeric(),
                  rmse=rmsez,
                  aic=aicz%>%as.numeric()%>%round(1),
                  n=nz%>%as.integer()
                )
              )
              
    }
    if(length(fez)==0){
      out_ <- data.table::data.table(
              model=modz[[m0]]["method"][[1]],
              y=ynamz,
              x=xnamz[1],
              ix=ifelse(is.null(ix),NA_character_,ix),
              beta_est=coef_t[1]%>%as.numeric(),
              beta_se=coef_t[2]%>%as.numeric(),
              beta_z=coef_t[3]%>%as.numeric(),
              beta_p=coef_t[4]%>%as.numeric(),
              beta_ixl_est = beta_ixl_est,
              beta_ixl_se = beta_ixl_se,
              beta_ixl_z = beta_ixl_z,
              beta_ixl_p = beta_ixl_p,
              beta_ixu_est = beta_ixu_est,
              beta_ixu_se = beta_ixu_se,
              beta_ixu_z = beta_ixu_z,
              beta_ixu_p = beta_ixu_p,
              converge=conv_status,
              offset=off_status,
              iv1_f=iv1_f,
              iv1_fp=iv1_fp,
              iv1_wu=iv1_wu,
              iv1_wup=iv1_wup,
              adj_r2=ar2%>%as.numeric(),
              pseudo_r2=pr2%>%as.numeric(),
              rmse=rmsez,
              aic=aicz%>%as.numeric()%>%round(1),
              n=nz%>%as.integer()
            )
    }
    out_
  }) %>% data.table::rbindlist(fill=TRUE)
  return(out)
}

# Labels
labmat <- data.table::data.table(
  varz=c("ln_N_EVENT_R","ln_N_ARREST_R","ln_N_ARREVT_R","ln_N_EVENT","ln_N_ARREST","ln_N_ARREVT","LONG","LAT","POP_WR","POP_W","TROIKA_DIST","KREMLIN_DIST","DIST2RAIL","DIST2STN","DIST2MTRO","DIST2GOVE","DIST2INDU","DIST2MUSE","DIST2EDUC","DIST2HOTE","DIST2ENTE","DIST2HOSP","DIST2TRAN","DIST2AGRI","DIST2NKVD","DIST2PROS","DIST2COUR","DIST2PART","DIST2MILI","DIST2CARC","DIST2CAMP","DIST2SHAR","DIST2CEME","N_TOTAL_R","log(N_TOTAL_R + 1)","PLA_PCTR","DENREMREQ_PCTR","I(PLA_PCTR/100)","I(DENREMREQ_PCTR/100)","DENREMREQ_PCTR/100","N_TOTAL","log(N_TOTAL + 1)","PLA_PCT","DENREMREQ_PCT","I(PLA_PCT/100)","I(DENREMREQ_PCT/100)","DENREMREQ_PCT/100","log(N_EVENT_R + 1)","log(EVT_PCTR + 1)","log(N_EVENT + 1)","log(EVT_PCT + 1)","fit_log(N_EVENT_R + 1)","fit_log(EVT_PCTR + 1)","fit_log(N_EVENT + 1)","fit_log(EVT_PCT + 1)","RAY_NAME","BLOK_ZONE","model","y","x","beta_est","beta_se","beta_z","beta_p","iv1_f","iv1_fp","iv1_wu","iv1_wup","fe1","fe2","fe3","fe4","fe5","fe6","adj_r2","rmse","aic","n","feols","feglm","rdd","sar","car","sp_de","sp_ie","sp_range","sp_extra","pseudo_r2","RESI_ID","BLOK_ID","nationality_short","party","okonkh","log(LA_ANY + 1)","LA_ANY","log(LA_DENREMREQ + 1)","LA_DENREMREQ","log(LA_ANY_2 + 1)","LA_ANY_2","log(LA_DENREMREQ_2 + 1)","LA_DENREMREQ_2","converge","offset","LPOP","INDU_DIST","AREA","log(N_PLAQO_Military_R + 1)","log(N_PLAQO_Party_R + 1)","log(N_PLAQO_Government_R + 1)","log(N_PLAQO_Arts_R + 1)","log(N_PLAQO_Science_R + 1)","log(N_PLAQO_Medicine_R + 1)","log(N_PLAQO_Exploration_R + 1)","log(N_PLAQO_Sports_R + 1)","log(N_PLAQO_Repressed_R + 1)","log(N_PLAQO_Other_R + 1)","log(N_PLAQO_Military_Sov_R + 1)","log(N_PLAQO_Party_Sov_R + 1)","log(N_PLAQO_Government_Sov_R + 1)","log(N_PLAQO_Arts_Sov_R + 1)","log(N_PLAQO_Science_Sov_R + 1)","log(N_PLAQO_Medicine_Sov_R + 1)","log(N_PLAQO_Exploration_Sov_R + 1)","log(N_PLAQO_Sports_Sov_R + 1)","log(N_PLAQO_Sov_Repressed_Sov_R + 1)","log(N_PLAQO_Other_Sov_R + 1)","log(N_PLAQO_Lenin_Sov_R + 1)","log(N_PLAQO_PartyNoLenin_Sov_R + 1)","log(N_PLAQO_Writer_Sov_R + 1)","log(N_PLAQO_Finearts_Sov_R + 1)","log(N_PLAQO_Architect_Sov_R + 1)","log(N_PLAQO_Music_Sov_R + 1)","log(N_PLAQO_Film_Sov_R + 1)","log(N_PLAQO_Film_R + 1)","log(N_PLAQO_Lenin_R + 1)","log(N_PLAQO_PartyNoLenin_R + 1)","log(N_PLAQO_Writer_R + 1)","log(N_PLAQO_Finearts_R + 1)","log(N_PLAQO_Architect_R + 1)","log(N_PLAQO_Music_R + 1)","nfirms_b","employees_b","corpsize_b","foreignowned_b","boardind_b"),
  labz=c("Executions","Executions","Event Size","Executions","Executions","Event Size","Longitude","Latitude","Population","Population","Distance to Troika","Distance to Kremlin","Distance to railroad","Distance to rail station","Distance to metro","Distance to government facility","Distance to industrial site","Distance to museum","Distance to educational facility","Distance to hotel","Distance to entertainment venue","Distance to hospital","Distance to transporation","Distance to agriculture","Distance to NKVD","Distance to prosecutor","Distance to court","Distance to party organ","Distance to military base","Distance to carceral facility","Distance to camp","Distance to sharashki","Distance to cemetery","Petitions","Petitions","Plaques (\\%)","Denied/Removed (\\%)","Plaques (\\%)","Denied/Removed (\\%)","Denied/Removed (\\%)","Petitions","Petitions","Plaques (\\%)","Denied/Removed (\\%)","Plaques (\\%)","Denied/Removed (\\%)","Denied/Removed (\\%)","Executions (logged)","Executions per capita (logged)","Executions (logged)","Executions per capita (logged)","Executions (logged)","Executions per capita (logged)","Executions (logged)","Executions per capita (logged)","Rayon","Zoning","Model","Outcome","Treatment","Estimate","Std. error","Z ccore","p value","First stage $\\mathcal{F}$ score","$\\mathcal{F}$ p value","Wu-Hausman","Wu-Hausman p value","FE 1","FE 2","FE 3","FE 4","FE 5","FE 6","Adj. R$^2$","RMSE","AIC","N","Linear","Binomial","FRDD","SAR","CAR","$\\sigma^2_{de}$","$\\sigma^2_{ie}$","$\\hat{\\rho}$","Extra parameter","Pseudo R$^2$","City block","City block","Ethnicity","Party","Industry","Petition","Petition","Denied/Removed","Denied/Removed","Petition","Petition","Denied/Removed","Denied/Removed","Converged","Offset","Population (log)","Distance to industrial site","Area","Military","Party","Government","Arts","Science","Medicine","Exploration","Sports","Repressed","Other","Military","Party","Government","Arts","Science","Medicine","Exploration","Sports","Repressed","Other","Lenin","No Lenin","Writers","Fine arts","Architects","Musicians","Film","Film","Lenin","No Lenin","Writers","Fine arts","Architects","Musicians","Number of tenants","Number of employees","Corporation size","Foreign-owned","Board independence"))
labmat_cat <- data.table::data.table(
  varz=c("industry","agriculture","forestry","transcomms","construction","retailfood","logistics","procurement","infoaccount","realestate","commerce","geodesic","othermaterial","residential","service","health","education","cultureart","science","econfinance","government","socialorganizations","foreign","unemployed","ARM","BEL","CHE","CHI","EST","GER","GRE","JEW","KAB","KAL","KOR","LAT","LIT","OSS","POL","RUS","TAT","UKR","higher","secondary","primary","illiterate","unknown","other"),
  labz=c("Heavy industry","Agriculture","Forestry","Transport/Communications","Construction","Retail/Food","Logistics","Procurement","Information/Accounting","Real Estate","Commerce","Geodesy","Other","Residential","Service","Health","Education","Culture/Arts","Science","Econ/Finance","Government","Social Organizations","Foreign","Unemployed","Armenian","Belarusian","Chechen","Chinese","Estonian","German","Greek","Jewish","Kabardin","Kalmyk","Korean","Latvian","Lithianian","Ossetian","Polish","Russian","Tatar","Ukrainian","Higher","Secondary","Primary","Illiterate","Unknown","Unknown")
  )

# LaTeX tables
xtabz <- function(modmat,yz,xz,modz,yordz,xordz,capz,labbz,fit_width=TRUE,subcaption=" ",fn=tempfile(),drop_stats=NULL){
  empty_columns <- colSums(is.na(modmat) | modmat == "") == nrow(modmat)
  if(grepl("ols|glm",modz)){
    modmat_ <- modmat %>% data.table::copy() %>% .[grepl(yz,y)&grepl(xz,x)&grepl(modz,model),] %>% .[,y:=factor(y,levels=yordz)] %>% .[,x:=factor(x,levels=xordz)] %>% .[order(y,x)] %>% .[,x := x%>%match(labmat[,varz]) %>% labmat[.,labz]] %>% .[,y := y%>%match(labmat[,varz]) %>% labmat[.,labz]]  %>% .[,eval(grep("^fe",names(.))):=lapply(.SD,function(x){x%>%match(labmat[,varz]) %>% labmat[.,labz]}),.SDcols=grep("^fe",names(.))] %>% .[,model := model%>%match(labmat[,varz]) %>% labmat[.,labz] %>% paste0(1:.N,". ",.)] %>% .[,beta_est := paste0(smart_round(beta_est,2))]%>% .[,beta_se := paste0("(",smart_round(beta_se,2),")",{beta_p%>%findInterval(c(0,.001,.01,.05,.1,1))%>%c("***","**","*","$^{\\dagger}$","")[.]})] %>% .[,adj_r2:=round(adj_r2,2)] %>% .[,pseudo_r2:=round(pseudo_r2,2)] %>% .[,rmse:=round(rmse,2)]  %>% .[,.SD,.SDcols=which(!empty_columns)]  %>% .[,c("beta_z","beta_p","x","converge","offset",drop_stats):=NULL] %>% data.table::transpose(keep.names="rn") %>% .[,rn:=rn%>%match(labmat[,varz])%>%labmat[.,labz]] %>% .[grep("FE",rn),rn:=paste0(V1," FE")] %>% .[grep("FE",rn),eval(names(.)[-1]):=apply(.SD,2,function(x){x[x!=""&!is.na(x)]<-"\\checkmark";return(x)})%>%data.table::as.data.table(),.SDcols=names(.)[-1]] %>% data.table::setnames(as.character(.[2,])) %>% .[-2]
  lpos <- list(
    pos=list(grep("Estimate",modmat_[,1][[1]])-1,grep("FE$",modmat_[,1][[1]])[1]-1),
    command=c(paste0("\\cline{2-",ncol(modmat_),"} \\\\[-1.75ex]\n "),paste0("\\cline{2-",ncol(modmat_),"} \\\\[-1.75ex]\n"))
    )
  }
  if(grepl("sar|car",modz)){
    modmat_ <- modmat %>% data.table::copy() %>% .[grepl(yz,y)&grepl(xz,x)&grepl(modz,model),] %>% .[,y:=factor(y,levels=yordz)] %>% .[,x:=factor(x,levels=xordz)] %>% .[order(y,x)] %>% .[,x := x%>%match(labmat[,varz]) %>% labmat[.,labz]] %>% .[,y := y%>%match(labmat[,varz]) %>% labmat[.,labz]]  %>% .[,eval(grep("^fe",names(.))):=lapply(.SD,function(x){x%>%match(labmat[,varz]) %>% labmat[.,labz]}),.SDcols=grep("^fe",names(.))] %>% .[,model := model%>%match(labmat[,varz]) %>% labmat[.,labz] %>% paste0(1:.N,". ",.)] %>% .[,beta_est := paste0(smart_round(beta_est,2))]%>% .[,beta_se := paste0("(",smart_round(beta_se,2),")",{beta_p%>%findInterval(c(0,.001,.01,.05,.1,1))%>%c("***","**","*","$^{\\dagger}$","")[.]})] %>% .[,pseudo_r2:=round(pseudo_r2,2)] %>% .[,sp_range:=round(sp_range,2)] %>% .[,rmse:=round(rmse,2)] %>% .[,.SD,.SDcols=which(!empty_columns)] %>% .[,c("sp_de","sp_ie","sp_extra","beta_z","beta_p","x"):=NULL]  %>% data.table::transpose(keep.names="rn") %>% .[,rn:=rn%>%match(labmat[,varz])%>%labmat[.,labz]] %>% .[grep("FE",rn),rn:=paste0(V1," FE")] %>% .[grep("FE",rn),eval(names(.)[-1]):=apply(.SD,2,function(x){x[x!=""&!is.na(x)]<-"\\checkmark";return(x)})%>%data.table::as.data.table(),.SDcols=names(.)[-1]] %>% data.table::setnames(as.character(.[2,])) %>% .[-2]
  lpos <- list(
    pos=list(grep("Estimate",modmat_[,1][[1]])-1,grep("FE$",modmat_[,1][[1]])[1]-1),
    command=c(paste0("\\cline{2-",ncol(modmat_),"} \\\\[-1.75ex]\n "),paste0("\\cline{2-",ncol(modmat_),"} \\\\[-1.75ex]\n"))
    )
  }
  if(grepl("rdd",modz)){
    modmat_ <- modmat %>% data.table::copy()%>% .[grepl(yz,y)&grepl(xz,x)&grepl(modz,model),] %>% .[,y:=factor(y,levels=yordz)] %>% .[,x:=factor(x,levels=xordz)] %>% .[order(y,x)] %>% .[,x := x%>%match(labmat[,varz]) %>% labmat[.,labz]] %>% .[,y := y%>%match(labmat[,varz]) %>% labmat[.,labz]] %>% .[,model := model%>%match(labmat[,varz]) %>% labmat[.,labz] %>% paste0(1:.N,". ",.)] %>% .[,beta_est := paste0(round(beta_est,2))]%>% .[,beta_se := paste0("(",round(beta_se,2),")",{beta_p%>%findInterval(c(0,.001,.01,.05,.1,1))%>%c("***","**","*","$^{\\dagger}$","")[.]})]  %>% .[,adj_r2:=round(adj_r2,2)%>%abs()] %>% .[,pseudo_r2:=round(pseudo_r2,2)%>%abs()] %>% .[,rmse:=round(rmse,2)] %>% .[,iv1_f := paste0("",round(iv1_f,2),"",{iv1_fp%>%findInterval(c(0,.001,.01,.05,.1,1))%>%c("***","**","*","$^{\\dagger}$","")[.]})] %>% .[,iv1_wu := paste0("",round(iv1_wu,2),"",{iv1_wup%>%findInterval(c(0,.001,.01,.05,.1,1))%>%c("***","**","*","$^{\\dagger}$","")[.]})]  %>% .[,.SD,.SDcols=which(!empty_columns)]  %>% .[,c("beta_z","beta_p","iv1_wup","iv1_fp","x","converge","offset"):=NULL] %>% data.table::transpose(keep.names="rn") %>% .[,rn:=rn%>%match(labmat[,varz])%>%labmat[.,labz]] %>% .[grep("FE",rn),rn:=paste0(V1," FE")] %>% .[grep("FE",rn),eval(names(.)[-1]):=apply(.SD,2,function(x){x[x!=""&!is.na(x)]<-"\\checkmark";return(x)})%>%data.table::as.data.table(),.SDcols=names(.)[-1]] %>% data.table::setnames(as.character(.[2,])) %>% .[-2]
    lpos <- list(
      pos=list(grep("Estimate",modmat_[,1][[1]])-1,grep("FE$|^F |mathcal",modmat_[,1][[1]])[1]-1),
      command=c(paste0("\\cline{2-",ncol(modmat_),"} \\\\[-1.75ex]\n "),paste0("\\cline{2-",ncol(modmat_),"} \\\\[-1.75ex]\n"))
    )
  }
  xmat <- modmat_%>%xtable::xtable(caption=capz,label=labbz,align = paste0("ll",paste0(rep("c",each=ncol(.)-1),collapse=""))) %>% data.table::setnames(paste0("\\multicolumn{1}{c}{",names(.),"}"))
  for(r0 in c(1:nrow(xmat))){xmat[r0,] <- paste0("\\multicolumn{1}{c}{",xmat[r0,],"}")}
  xmat <- xmat%>%xtable::print.xtable(include.rownames=FALSE,add.to.row=lpos,booktabs=TRUE,sanitize.text.function=identity) %>%
      gsub("\\*\\*\\*","\\*\\*",.) %>%
      gsub(" 0 "," 0.0000 ",.,fixed = T) %>%
      gsub(" (0)"," (0.0000)",.,fixed = T) %>%
      gsub("(0)","(0.0000)",.,fixed = T) %>%
      gsub("\\multicolumn{1}{c}{NA}","\\multicolumn{1}{c}{}",.,fixed = T) %>%
      gsub("\\multicolumn{1}{c}{Denied/Removed (\\%)} & \\multicolumn{1}{c}{Denied/Removed (\\%)}","\\multicolumn{2}{c}{Denied/Removed (\\%)}",.,fixed = T) %>%
      gsub("\\multicolumn{1}{c}{Petitions} & \\multicolumn{1}{c}{Petitions}","\\multicolumn{2}{c}{Petitions}",.,fixed = T) %>%
      gsub("\\multicolumn{1}{c}{Denied/Removed} & \\multicolumn{1}{c}{Denied/Removed}","\\multicolumn{2}{c}{Denied/Removed}",.,fixed = T) %>%
      gsub("\\multicolumn{1}{c}{Petition} & \\multicolumn{1}{c}{Petition}","\\multicolumn{2}{c}{Petition}",.,fixed = T)
  if(fit_width==TRUE){
    xmat <- xmat %>% gsub("\\begin{tabular}","\\makebox[.95\\textwidth]{\\begin{tabular}",., fixed = TRUE) %>% gsub("\\end{tabular}","\\end{tabular}}",., fixed = TRUE)
    out <- cat(gsub("\\end{tabular}}\n", paste("\\end{tabular}}\n\\caption*{{\\footnotesize ", subcaption, "\n}}\n", sep = ""), xmat, fixed = TRUE), sep = "\n", file =  fn)
  }
  if(fit_width==FALSE){
    out <- cat(gsub("\\end{tabular}\n", paste("\\end{tabular}\n\\caption*{{\\footnotesize ", subcaption, "\n}}\n", sep = ""), xmat, fixed = TRUE), sep = "\n", file =  fn)
  }
  return(xmat)
}

# BMA table LaTeX
xtabz_bma <- function(bma_mod_1,bma_mod_2,labmat,capz,labbz,festart,fit_width=TRUE,subcaption=paste0("$P(\\hat{\\beta}\\neq0)$ are posterior inclusion probabilities, scaled 0 to 100. Quantities reported under ``Post. Mean'' are means (standard deviations) of coefficient estimates' posterior distributions. ",bma_mod_1["n.models"][[1]], " Gaussian models, ",bma_mod_2["n.models"][[1]], " quasi-Poisson models selected (",bma_mod_2["n.vars"][[1]]," covariates)."),fn=tempfile()){  
  options(scipen=999)
  bma_mat <- data.table::data.table(
    Variable=c("Intercept",bma_mod_1["namesx"][[1]]),
    probne0_1=c(100,bma_mod_1["probne0"][[1]]),
    postmeansd_1=bma_mod_1["postmean"][[1]] %>% smart_round(1) %>% paste0(" (",bma_mod_1["postsd"][[1]] %>% smart_round(1),")"),
    probne0_2=c(100,bma_mod_2["probne0"][[1]]),
    postmeansd_2=bma_mod_2["postmean"][[1]] %>% smart_round(1) %>% paste0(" (",bma_mod_2["postsd"][[1]] %>% smart_round(1),")")
    ) %>% .[Variable%in%labmat[,varz],Variable:=Variable%>%match(labmat[,varz])%>%labmat[.,labz]] %>% data.table::setnames(c("Variable","$P(\\hat{\\beta}\\neq0)$ (1)","Post. Mean (1)","$P(\\hat{\\beta}\\neq0)$ (2)","Post. Mean (2)"))
  lpos <- list(
      pos=list(grep("Intercept",bma_mat[,1][[1]])-1,grep(festart,bma_mat[,1][[1]])[1]),
      command=c(paste0("\\cline{2-",ncol(bma_mat),"} \\\\[-1.75ex]\n "),paste0("\\cline{2-",ncol(bma_mat),"} \\\\[-1.75ex]\n"))
    )
  xmat <- bma_mat%>%xtable::xtable(caption=capz,label=labbz,align = paste0("ll",paste0(rep("c",each=ncol(.)-1),collapse=""))) %>% data.table::setnames(paste0("\\multicolumn{1}{c}{",names(.),"}"))
  for(r0 in c(1:nrow(xmat))){xmat[r0,2:ncol(xmat)] <- paste0("\\multicolumn{1}{c}{",xmat[r0,2:ncol(xmat)],"}")}
  xmat <- xmat%>%xtable::print.xtable(include.rownames=FALSE,add.to.row=lpos,booktabs=TRUE,sanitize.text.function=identity) %>%
      gsub("\\multicolumn\\{1\\}\\{c\\}\\{Variable\\}","\\multicolumn{1}{l}{Variable}",.)%>%
      gsub("\\*\\*\\*","\\*\\*",.) %>%
      gsub("\\$ \\(1\\)","$",.) %>%
      gsub("\\$ \\(2\\)","$",.) %>% 
      gsub("Mean \\(1\\)","Mean",.) %>%
      gsub("Mean \\(2\\)","Mean",.) %>% 
      gsub(" 0 "," 0.0000 ",.,fixed = T) %>%
      gsub(" (0)"," (0.0000)",.,fixed = T) %>%
      gsub("(0)","(0.0000)",.,fixed = T) %>%
      gsub("\\multicolumn{1}{c}{NA}","\\multicolumn{1}{c}{}",.,fixed = T) %>%
      gsub("\\toprule\n","\\toprule\n&\\multicolumn{2}{c}{(1) Gaussian}&\\multicolumn{2}{c}{(2) Quasi-Poisson}\\\\",.,fixed = T)
    if(fit_width==TRUE){
      xmat <- xmat %>% gsub("\\begin{tabular}","\\makebox[\\textwidth]{\\small\\begin{tabular}",., fixed = TRUE) %>% gsub("\\end{tabular}","\\end{tabular}}",., fixed = TRUE)
      out <- cat(gsub("\\end{tabular}}\n", paste("\\end{tabular}}\n\\caption*{{\\footnotesize ", subcaption, "\n}}\n", sep = ""), xmat, fixed = TRUE), sep = "\n", file =  fn)
    }
    if(fit_width==FALSE){
      out <- cat(gsub("\\end{tabular}\n", paste("\\end{tabular}\n\\caption*{{\\footnotesize ", subcaption, "\n}}\n", sep = ""), xmat, fixed = TRUE), sep = "\n", file =  fn)
    }
    options(scipen=5)
    return(xmat)
}

